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ABSTRACT 


Several areas of applications of image processing to astronomy have 
been addressed in this report. These areas include: 1) polar coordinate 

transformations: a means of enhancing the faint outer spiral pattern of 

M83; 2) ra!:los of spectral bands in M82 indicate that the Ha filaments 

are not of a reflection nature; 3) extending the Gull maximum entropy 
algorithm beyond the convergence point produces better restoration with 
little artifact. Applications to M87; 4) attempts to recover Images 

from the phase of the Fourier transform produce uncertain rasults in 
real objects. Applications to M87; 5) a proposed algorithm to statis- 

tically identify more unique morphological types of galaxies is illustrated 
based upon the correlation of images; 6) modeling of complex nebulosity 
In Lick Ha 101 with a median filter in order to automatically catalogue 
stars with a threshold algorithm; 7) relative color and polarization In 
R-Monocerotls Indicate a hollow, cone shapeJ, dusty reflection nebula; 

8) CCD type aliasing. Periodic scenes can be corrupted by aliasing if 
the image forming optics are of too high a quality; 9) illusions induced 
by sampling frequency beating. These illusions can be eliminated by 
doubling the picture size with the Sampling Theorem, and 10) low dynamic 
range induced illusions. Unexpected structure can appear in an FFT If 
the dynamic range of the image Is small. 
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SECTION I 


INTRODUCTION 


Astronomy has traditionally made great use of Imagery, but this data 
has almost exclusively been evaluated via the photographic process. With 
the ndVv'Ot of the unmanned planetary missions and digital cameras, tech- 
niques for processing digital Images have been developed and are beginning 
to be applied in astronomy. The advantages of digital images over photo- 
graphic ones are very great and, with the development of new algorithms and 
cameras, become greater all the time. It is thus of importance to search 
for applications and disseminate digital techniques to the astronomical 
community at large. When Space Telescope becomes operational, the quantity 
and unprecedented quality of the calibrated digital images returned will 
call for every available processing means at our disposal. This report is 
Intended to prov-ide a few of those means. 


This Work was performed during 1980 at the Image Processing Labora- 
tory (IPl.) of JPL under FY80 RTOP f/389-4l-0l. It represents a continuing 
effort (References 1, 2, and 1) to Identify and develop now image process- 
ing applications to astronomy. Additional goals of this research are the 
further education of the astronomicii I community in digital image process- 
ing techniques and the expansion of methods for analysis of high-quality 
space telescope imagery when it becomes available. Hopefully, others will 
be inspired to find applications to their own disciplines. 


Tliis report is divided into twelve sections. Sections II 
present techniques applied to a major topic. The data used and 
nlques involved are described in detail, along with photographs 
of the results. Section XII gives conclusions that can be made 
of this work. 


through XI 
the tech- 
and plots 
as a result 
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SECTION II 


POLAR COORDINATE TRANSFORMATIONS - M83 


It is often ndvantageouii to tranti'orm Images geometrically Into 
other than Cartesian coordinate systems. Two reasons are commonly given 
for this (T’.eference 4): 

1. To take advantage of symmetry In the object. 

2. To generalize the effect of an otherwise Cartesian algorithm. 

As an Illustration of such a process, deep plates of NS) obtained 
by Dr. Willem Wamsteker of the ESO were selected. Figure 1 Illustrates 
M63 and what appear to be spiral arms In a ring surrounding the object. 

In order to reveal faint azimuthal surface brightness It would be desir- 
able to smear the object In rotation, Improving the signal to noise for 
low frequency azimuthal structure. Conventional filters are usually coded 
to operate In a Cartesian system only and would not normally be of use to 
us In this case. If, however, we vere to transform the image to polar 
coordinates such that each line represented a different radical profile, 
extending In azimuth from 0" to 360° we would be able to apply a rotational 
filter by merely o|>eratlng upon the polar coordinate Image. A subsequent 
inverse transformation would return us to the original gec.netrl • 
representation. 

Figure 2 displays the polar coordinate Image of M83 performed after 
a 41 by 41 element median filter (Reference 4) to suppress small stars. 

The left edge of the picture Is the galaxv nucleus, the lines extending 
from the top down clockwise about the nucleus from the right. Figure 3 
Illustrates N83 after retransforming to C.irteslan coordinates again. In 
Figure 3, M83 has been rotatioaally smeared by 30°. This was acccmipllshed 
by operating a vertical median filter of length 30/360 NL on the polar 
coordinate image, where NL is the number of lines in the image. Several 
large stars were not eliminated by the filter and have been smeared, A 
faint outer spiral arm to the upper right of M83 can now be followed con- 
siderably further than In the original image. 
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Figure 1. M83 contra*<t enhanced to dlaplay the background aky and 

large-scale emulalon Inhonuigenettlea. Extended faint 
azimuthal spiral arms are visible. 
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Figure 2 , M83 polar coordinate tranti* 

formation. The lines corre- 
spond from top to bottom with 
radial profiles clockwise 
about the galaxy nucleus 
beginning from the right. 

The nucleus lies along the 
left edge. A 40 by 40 
median filter was applied 
before the transformation. 
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Figure 3. N83 after retransformlng from polar coordinates and 

after rotational smearing of 30 degrees with a median 
filter to Improve the azimuthal signal to noise at low 
spatial frequencies. 
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SECTION III 


MULTISPECTRAL RATIOS - M82 


Plates of M82 obtained by Dr. Alan Sandaga and procaaaad pravioualy 
for color and polarization (Rafaranca 3) have bean ratioad to dataraina 
tha natura of tha Ha aoiirca. Tha two thaoriaa for tha Ha aourca era 
raflaction of Ha from tha nuclaua off of tha filaaanta (Rafaranca 5), and 
emiaaion of Ha from tha filamanta .thamaalvaa (Rafaronca 6). Figuraa A, 5, 
and 6 llluacrata tha Ha/broad band ratios in the aanae Ha/103a0 (blua)» 
Ha/103aD (graan) , and Ha/103aE (rad) raapactivaly. From an inapactlon of 
Piguraa 4-6 it ia apparent that tha ratio of Ha/broad band ia a atrong 
function of color, being tha moat dlatinct in tha blue. Thla would argue 
againat tha aourea of H being raflaction off of duat in the filamanta. 

If tha aourca of all filament Irradianca ware tha nucleus than tha fila- 
ments should remain constant in all tha ratios. 

The white regions to upper right and lower left in Figures 4-6 are 
due to division by small numbers close to the sky background and are not 
real. 
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Figure 4. M82 ratio Image of Ha/103a0. 
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Figure 5. M82 ratio image of Ha/103aD . 
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Figure 6. M82 ratio image of Ha/103aE. 
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SECTION IV 


EXTENDED MAXIMUM ENTROPY RESTORATION - M87 


In Reference 3» results were presented for the maximum entropy 
restoration algorithm developed by Dr, Stephen Gull (Reference 7). It 
was shown that the algorithm possessed the valuable quality that It Intro- 
duced a minimal amount of ringing; however. It also was so conservative 
as to produce only a slight improvement in resolution. The algorithm 
has been modified to continue the restoration process beyord the normal 
convergence point by Incrementing the qu.’.ntlty \ by successive multiplica- 
tion. In Figure 7 a noiseless target (upper left) was blurred by an 
11 by 11 element box filter and restored by the extended maximum entropy 
technique (center left) and, for comparison, by a Wiener filter with a 
signal to noise of 15 (center right). The restorations are now virtually 
Identical, a considerable improvement over the old algorithm dlscusred 
In Reference 3. To the lower left Is the chl-squared residuals map with 
the MTF of the 11 by 11 filter to Its right. 

In Figure 8 a test target has been blurred (upper right) and then 
decimated (lower left). The extended maximum entropy algorithm has been 
able to restore the characters to a legible level by ignoring missing 
pixels. It would not be possible to operate a linear filter on such an 
image . 


For application to a real object, images of M87 obtained by Drs. Arp 
and Sulentic were used. In Figure 9 each row contains the M87 original 
image to the left, the extended maximum entropy restoration made from the 
image to its left at center, and the chl-squared residuals to the right. 
The M87 iniagu! In the left lower row is the result of subtracting a model 
of the skv + galaxy halo from it in order to aid the non-negativity con- 
straint. considerable resolution gain of about 80% is evident, most 
noticeable in the lover row. 

The restored images suffer from little artifact and are superior to 
linear restorations. This algorithm should be preferred to linear filters 
when the true nature of the object is unknown. 
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Figure 7. Restoration of a noiseless blurred target; (upper left) 
original; (upper right) blurred with an 11 by 11 bi>x- 
f 11 ter; (center left) extended maximum entropy; (center 
right) Wiener; (lower left) maximum entropy chl-squared 
residuals; (lower right) transfer function (MTF) of the 
boxfllter. The Wiener signal to noise was 15. 
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Figure 8. Restoration of a blurred and decimated Image; (top left) 
original: (top right) blurred with a 7 by 7 boxf liter; 
(lower left) decimated revision of the top right Image; 
(lower right) extended maximum entropy restoration. 
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Figure 9. M87 moKalc. Each row displays the original linap,e at the left of 
the extended maximum entropy image in the center, and the chi- 
squared residuals to the right. A model of the sky and the M87 
halo was subtracted from the original image in the bottom row. 
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SECTION V 

RESTORATIONS FROM THE PHASE - M87 


It Is posalble to recover the amplitude from the phone of a Fourier 
transform and vice versa provided that the Inverse transform has certain 
properties (such as positivity; all of the necessary conditions are not 
known). A flow chart is shown h low which Illustrates the basic concept 
originally suggested by Flenup (Reference 8). 



The strategy is beginning with an initial random noise estimate of 
an object, to force either the amplitude or the phase of the FFT to match 
the amplitude or the phase of the true FFT (whichever is known). The 
inverse FFT of the estimate is corrected for negative intenSi-tles and the 
FFT is once again forced to match the reference FFT. The process continues 
until a satisfactory image is obtained (if that proves possible). The 
YES/NO test in the last step is used to "shock" the solution away from 
degenerating into a repetitive static situation. 

Figure 10 Illustrates this process using a test target (top center) 
and a typlca’ initial estimate for the target (top left). All of the 
other squares represent images recovered from only the phase of the FFT 
of the test target after 10 iterations, each with a different initial 
random noise estimate. These images are very similar, and retain the 
point spread function of the original target (note that any symmetrical, 
non-axis crossing modulation transfer function would have served). 


I 
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ThlB concept was applied to an image of M87 in Figure LI. The left 
and right columns contain the original "test" images at the top and the 
images recovered from the phase of their FFT below them. The left column 
used the original image and the right column used the original maximum 
entropy algorithm restoration. It was not expected that the loft column 
would produce an improved image because the phase of the FFT has been 
corrupted by film grain noise. In the right image however, the phase has 
been largely restored even though the amplitude has not been resolution 
corrected, making this a good candidate for this technique. The results 
show no improvement in the left column and a questionable result in the 
right column. 

The M87 Image may hove greatly improved in resolution, however, the 
row of reference stars below has been obliterated, a result which obviates 
any physical interpretation of the M87 image itself. The likely answer 
for failure in this case is initial errors in the phase. 

A possible correct'ive avenue of approacit might be to "rig" the image 
wit dummy structure convolved with the blurred image point spread function 
and then each time the inverse FFT is taken, to replace the dummy structure 
with the actual unconvolvcd structure in order to guide the restoration 
process in the correct direction. 
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Figure 10. 



Test target for restorations from the phase; (top left) 
typical random noise Initial estimate of the true scene; 
(top center) the true scene; all other Images are the 
product of 10 Iterations using the phase of the FFT 
only to retrieve the amplitude. 
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Figure 11. M87 mosaic of restorations from the phase; (top left) original M87 image; 

(top right) maximum entropy restoration without sky model subtraction; 
(lower left) restoration from the phase using the upper left image; 

(lower right) restoration from the phase using the upper right image. 
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SECTION Vi 


PRINCIPAL COMPONENT MORPHOLOGICAL TYPES - A PROPOSAL 


The following formAlizetion wee developed efter e diecueelon with 
Dr. Merlaue Molea Vlllemete of Centre Aetronomlco Illepeno-Aleiiien. 

The object of the following topic will be to eelect from e lerge 
eelection of galexiee repreeenting ell known oiorphologicel typee the 
etetieticelly fundementel morphologicel typee, i.e. , thoee which correlete 
the leeet with the ochere end ere thue the moet unique. In order to eccom- 
plieh thie« we muet compute the correletion metrix which coneiete of the 
correlation of each galaxy with all of the othera. 


•ij 




S Kj ■ ‘j) ’ (‘kl ■ u) ' ‘ 


< i < N 
1 < J < N 


( 1 ) 


where 

I > intensity of pixel k in image 1 or j 
1 ■ mean intensity of image i or J 
N ■ total number of images 
0 - standard deviation of image i or j 
C^j is an N by N matrix 

The eigenvalues of must be computed next. 

Once the el^^^nvalues have been computed the eigenvectors can be 
computed from: 

R*^CR - C' (2) 


where 

t ■ transpose 
R ■ eigenvector matrix 

C' ■ matrix whose diagonal is the eigenvalues 
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The principal componant tranaformation (Rafarancaa 3 and 9) of 
convantional imagaa into uncorralatad Imagaa la glvan by multiplying 
matrix R with tha convantional Imagaa: 


hi\ /‘i\ 

S " ** 

a a 

VI VI 


( 3 ) 


The elgenvaluea represent the information content (variance) of the 
uncorrelated imagaa I'. 

Thus, if there are I large eigenvalues and n-H smaller ones It means 
that only it unique Images exist In the original set of N Images !• and that 
the other n>il images arc merely linear combinations of the former l Images. 
This tells us how many unique morphological types exist but it does not 
specify which ones. Notice that for arbitrary pictures the unique Images 
may prove to be linear combination of the Image provided; however, we know 
from the physical positivity constrain on Intensities that virtually no 
linear combinations are possible which provide only positive Images. This 
means that the set of unique morphological types must be single images 
contained In the set provided. 

The unique Images can be Identified from the top row of the R matrix 
which is that vector of weights which, when multiplied by the Input set of 
galaxies, provides that image with the largest eigenvalue. The principal 
component image with the largest eigenvalue Is the optimal average of the 
Input Images, and consists of the weighted rows of the Input Images. The 
weighting Is large If an Image Is "typical", l.e. , Is highly correlated 
with Its neighbors, and Is small If It Is different, l.e., is unique among 
the set of input Images (a distinct morphological type). 

Figure 12 illustrates the process described thus far. The top row 
consists of four Images, only two of which are unique?. The third Image 
from the left Is the average of Images 1 and 2 and the fourth Image con- 
sists of 20% of image 1 plus 80% of image 2. The correlation matrix C . . 

Is: ^ 
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0.46 

0.85 

0.62 
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0.85 
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0.85 

0.85 
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0.94 

0.62 

0.98 
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Proffl C the elgenvaliieH can be computed aa: 


3.381 

0.620 

0.001 

O.OOl 


only the top two of which are significant, indicating that only two unique 
Images exist. 

The eigenvector matrix H ran be computed as: 



0.4282 

0.4948 

0.5401 

0.5293 \ 

-0.7831 

0.5277 

-0.1494 

0.2929 ’ 

0.268/ 

-0.1563 

-0.6998 

0.6431 

-0.3619 

-0.6724 

0.4428 

0.4696 j 


iCach of 
bottom row of 


the four Images (the principal component Images) on the 
Figure 12 are obtained from the R matrix: 




\ 


n 

Hi* 


bottom 



top 


From an Inspection of the top t*nw of the R matrix it can be seen 
that the left two values are the sma ast. Indicating that the left two 
Images In the top row are the most unique of the set of 4 Images we began 
with. 


From a selection of a hundred Images of all morphological types each 
face-on, Intensity rotation and size normalized It should be possible 
cither to select a few parent types or at least to rank the Images In 
Importance, 
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Figure 12. Demonstration target Illustration of the principal 

component transformation. The upper four Images are 
linear combinations of two Images. The lower four 
Images are normalized principal component Images of the 
upper four Images ranked from left to right In decreas- 
ing eigenvalue. From the lower Images we can determine 
which of the upper four are the original two Images. 
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SECTION VII 


AUTOMATED COMPUTATION OF STELLAR MAGNITUDES IN 
NEBULOSITY - LICK Ha 101 


One of the problems encountered In automated methods of determining 
stellar magnitudes Is when spatial nebulosity causes the background of the 
plate to vary. This situation la a serious drawback for stellar catalog- 
ing software vhlch relies upon a simple above- threshold criterion for object 
detection. The median filter offers an expedient If not a very accurate 
solution to this problem by providing a means of modelling the background 
without the stars. Once the background Is determined It can be subtracted 
from the Image, leaving the stars on a low amplitude Intermediate spatial 
frequency background. 

Figure 13 Illustrates one of a set of multispectral Images of Lick 
Ha 101 after the background has been subtracted and the stars catalogued. 
Considerable confusion exists near the complicated portions of the nebula 
resulting In numerous false "stars” entered Into the catalogue. From the 
numerical designations It Is however quite easy to Identify reasonable 
objects and to select their Intensities from among the sorted catalogue 
entries. 
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Figure 13. Lick Ha 101 stellar catalogue map. Stellar images were 
located by simple threshold search after subtraction of 
a model of the Lick Ha 101 reflection nebula. The 
nebula was modelled with a median filter. 




SECTION VIII 


COLOR AND POLARIZATION IN R-MONOCEROTIS 


The enhancement of color and the generation of polarization maps Is 
not new and has been discussed In References 1 and 3. This section Is 
Included because, for R-Monocerotls, both are available together and the 
combination provides an Interesting data base on an early type object. 
Figure lA Is a relative color display produced by normalizing the blue, 
green, and red plates (103a0, 103aD, 098) and exaggerating the color with 
a saturation Image contrast enhancement xn the HSI space (Reference 4). 
Figures 15 and 16 Illustrate polarization (A percent and above has been 
saturated to white) and polarization vector orientation respectively. 

From the vector map It would appear that the large R-Monocerotls nebula 
Is a dusty reflection nebula Illuminated by the bright star at Its apex. 

The relative color display shows that the outer reflection nebula Is redder 
than at Its center, probably because the light Is more highly reddened 
by scattering at the periphery of the nebula. This would be the case If 
the nebula were a hollow cone. Light scattered at the sides of the cone 
would have to travel through more material before exiting toward the 
observer than would light at the front and back sides. The cone hy- 
pothesis Is further strengthened because the polarization map Is corre- 
lated with color, being more highly polarized where the nebula Is redder. 

Similar maps of the T-Taurl reflection nebula show no color varia- 
tions whatsoever. The different spectral bands are multiplicative replicas 
of each other. 
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Figure 14. R-Monocerotls reflection nebula. The three-color bands 
were normalized to the same amplitude to reveal differ- 
ential color. The true color is entirely red. 
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P'iguro 15. R-Monocerot Is polarization codod ft>r brightnoss. Four 
peroont polarization and above Is saturated. The 
vertical trend In the background is due to shading In 
one of the Images and is not real. 
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SECTION IX 


CCD ALIASING 


The following three sections were developed Jointly with Dr. Alan 
Gillespie of JPL and appear In more complete form as Reference 10. 

The sampling theorem (References 11 and 12) states that a band-limited 
function whose Fourier transform Is zero at frequencies greater than some 
frequency S^. Is fully specified by values spaced at equal Intervals not 
exceeding (2 with the exception of harmonic terms with zeroes at 

the sampling points. Real scenes have extended spectra» but may be band- 
limited by filtering In the atmospherei the optics* or the sensor Itself. 
CCD's and other solid-state sensors cannot band-limit a signal properly 
because the sample dimension Is at most equal to the sample spacing. Thus 
when the point-spread function (PSF) of the combined atmosphere and optics 
Is similar in size to Individual receptors of the CCD* care must be taken 
in interpreting resulting Images. 

The process of sampling Is briefly summarized In Figure 17. A con- 
tinuous scene might have an extensive spectrum F(s) as shown In Figure 17a. 

A CCD or other solid-state camera with a square receptor (pixel) of width 
d will have a transfer function specified by sin (Hds)/]Ids* or slnc(ds)* 
where s Is the spatial frequency. The energy spectrum available to the 
CCD for measurement Is (F(s) 'slnc (ds) * a sine function with attenuated side 
lobes shown In Figure 17b. When the scene Is sampled at an Interval p* 
the continuous spectrum Is converted to a repetitive and overlapping array 
of discrete spectra with a periodicity of p~l In the frequency domain. For 
CCD cameras* d £ p and the fundamental spectrum Is significantly contamin- 
ated by energy from adjacent spectra (Figure 17c). For example, at 
0.25 cycles/pixel* half the energy may be aliased. If d > 2p, the only 
contamination would be from side lobes of neighboring spectra. Figure 17d 
Illustrates such a spectrum* specified by: 


00 

1(3) * XI *'p) ■ [*' (® "^f)] 

n—00 


(4) 


Even so, as evaluation of the above series shows, roughly one quarter (the 
exact amount depends on F(s)) of the energy at 0.25 cycles/pixel Is aliased. 

In an Image, aliased energy Is manifest as fine detail which Is dis- 
torted In two ways: small features may appear too large and their loca- 
tions may be shifted by up to d//F. As Root (Reference 13) has Illustrated, 
for most scenes the distortions caused by CCD - style undersampling are not 
severe and Indeed may Improve the appearance of the Image by Increasing the 
"sharpness" and the amount of detail visible. This Is because of the large 
amount of energy aliased when d ■ p. Most such Images can be Interpreted 
correctly, provided the scene Is familiar to the analyst. Simple geometric 
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targets may be grossly distorted, however, and scenes containing repetitive 
geologic structure or periodic atmospheric cloud features may likewise be 
misrepresented. It should be emphasized that low-pass filtering can do 
little to reduce these artifacts If any scene detail Is to be preserved at 
all, because fundamental and aliased energy cannot be separated and are 
both reduced alike. Convernely, high-pass filtering emphasizes the most 
strongly aliased part of the spectrum, exacerbating interpretative 
difficulties. 

Figures 18 Illustrates the visual impact of CCD-type undersampling 
on n highly structured target. The original "scene" is at the upper left. 
The picture at the upper right represents the scene band-limited such that 
sampling at the CCD pixel spacing p produces Isolated periodic Image spec- 
tra, as In Figure 17d. For comparison, at center left Is the image result- 
ing from sampling the original "scene" with d ■■ p simulating the image 
from a CCD camera. At bottom are the log amplitude of the Fourier trans- 
forms of the top Images. Finally, at right center Is the CCD Image after 
post-sampling low-pass filtering, illustrating the Inf easlblllty of this 
approach as a remedy for undersampling. 

The essential point Is that in the presence of quasl-perlodlc 
structure, image information can be so fundamentally altered by CCD-type 
undersampling and aliasing as to preclude correct visual Interpretation. 

In the properly sampled image (upper right) Information in the overlapped 
spectral region is removed, protecting the analyst from false interpreta- 
tions by rendering none possible. Further examples of sampling artifacts 
with various d and p may be found in Reference 13. 

Once a scene has been undersampled there is no simple procedure for 
resampling or filtering the data to produce a non-allased image. This is 
because phase information has been Irrevocably lost during integration for 
each pixel. Thus such intuitively simple and appealing suggestions as 
addition of neighboring pixels to simulate original sampling with d ■ 2p 
will fai] , at least in detail. Addition of adjacent pixels is a convolu- 
tion operation specified In one dimension by the kernel (0 1 1), which has 
the Fourier transform 


G(s) 


St 


-ilTS 
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COSTTS 


(5) 


when d “ p “ 1. Applied to the image spectrum, this filter modifies the 
spectrum amplitude by cos ns and also shifts the phase by yielding 

a filtered image spectrum of 


00 

i TTs 

e cosns y F (s + n) sine (s + n) (6) 
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which l8 not at all the same as the Intended spectrum specified by 


ao 

F (s + n) sine [2(8 + n)] (7) 

na-o> 


This basically reflects the Impossibility of separately filtering the In- 
dividual of the series once summation has occurred during sampling. 

Aliasing In CCD images Is most readily solved by degrading the qual- 
ity of the Image-forming optics. This has the effect of reducing the 
amplitude of F(s) at high frequencies so that overlap of the spectra of 
the sampled scene is negligible. This requires, however, that the optical 
PSF be approximately gaussian in order to avoid Introduction of optical 
artifacts. Simple defocusslng Is not a solution because of the onset of 
spurious resolution In the center of the modulation transfer function at 
defocuslngs greater than about 5.6X f2, where X Is the wavelengtii of light 
and £ is the optics f number (Reference 14). For example, an optical sys- 
tem operating at f/lS and 5000 A would begin to suffer phase reversals of n 
at a defocus greater than 0.6 mm. 
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Figure 17. Summary of Hampllng proceBBea: (a) F(s) represents the 

continuous Fourier Transform of a scene projected onto 
the Image plane containing the sensor. Frequency Is 
given by 8. In this Instance, the "lens" was assigned 
an Impulse PSF; real optics might have a sine or gauss Ian 
PSF, resulting In decreased amplitude of F(s) at large s; 
(b) Integration of the image photo receptors of width d 
modifies the image spectrum; (c) After sampling the Inuige 
at intervals of p the discrete Image spectrum Is the sum 
of replicates of the continuous spectrum repeated at 
frequency Intervals of p"^ cycles/pixel. For CCD cameras, 
d < p and severe overlapping will occur; (d) If D « 2p or 
If the continuous image spectrum Is band- 11ml ted, over- 
laps of spectrum replicates can be minimized. 
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Igure 18. Effects of CCD-style undersampling on a simple sinusoidal 
spoke target. Upper right: properly sampled Image 

(d ■ 2p) of the band-limited target showing absence of 
Information and artifacts at high frequencies. Center 
left: target sampled as seen by CCD (d ■ p) . Artifacts 

caused by aliasing are seen Instead of high-frequency 
data. Center right: CCD-sampled Image after band-pass 

filtering showing Impossibility of artifact removal. 

Lower left: spectrum of target; lower right: spectrum 

of correctly sampled image. 
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SECTION X 


SAMPLING FREQUENCY BEATING INDUCED ILLUSIONS 


Even If a scene Is band- limited and properly sampled* there la no 
guarantee that It will be recognized from visual Inspection of an arbitrary 
pictorial representation of the digital data. The sampling theorem only 
states that exact reconstruction of the function Is possible. In order 
to do this correctly* It Is necessary to Interpolate using the Infinite 
sine function. This Is not easy to do* and usually sampled data Is dis- 
played without any effort to reconstruct the original function. 

One consequence of failure to reconstruct the continuous Image before 
display Is the visible presence* obvious In ordered Images* of Interference 
patterns resulting from sampling harmonic constituents of the scene near 
their nodes* where the amplitude Is zero. The frequencies of these nodal 
zones can be specified In one dimension by observing that the <:ampllng 
function If given by 


6 (x + n) ■ jcostrx| (8) 

Harmonic constituents of the spectrum are waves in space described by a 
cos 2tt(u)x + <})) , where A Is amplitude. ui Is frequency, x Is spatial posi- 
tion* and <p Is phase. At some Xq the wave will be sampled at its crest. 

If Xq Is taken as the origin* the phase at which sampling occurs Is seen 
to be the periodic function 


<()(x) - X (0.5 - w) 


(9) 


To the e>?, the characteristic frequency of (J)(x) in an image is set by the 
repetition of the nodal zones, where ■ ± n/2. These occur twice per 
wavelength, so the apparent frequency of the nodal patterns Is l-2u. For 
u) < 0.25 this nodal frequency exceeds the Nyquist (0.5) and Is aliased. 

Figure 19 illustrates five waves of different to. Each is oversam- 
pled and should resemble a simple sinusoid. Wavelengths from top to 
bottom are 40, 40/6, 40/10, 40/15, and 40/19 samples. A wave of the sam- 
pling frequency has a period of 40/20 samples and Is represented as a 
binary sequence whose amplitude may range frc' ±1 to ±0 as a function of 
phase. It is evident that the simple nature of the sinusoid Is masked for 
wavelengths shorter than about four samples, and that Interference pat- 
terns may actually be the prominent feature of the image. 

Visual inspection of images dominated by waves of frequencies greater 
than 0.25 cycles/pixel will result in Inaccurate interpretation, even though 
the sampling theorem was observed. Because a discrete frequency spectrum 
of each of these sampled waves is simply II (s) » 6 (s + o)) + 6 (s - oj) , 
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where 6 Is the Imptilse function, s Is frequency, and (o is the frequency of 
the wave. It follows that the nodal or Interference patterns are com- 
pletely Illusory. In fact, they result from incorrect display of the dis- 
crete Image. 

The proper way to reconstruct the continuous Image requires convol- 
ving the Image by sine x and resampling at infinitesimally small Intervals 
p. Clearly this cannot be done, but correct resampling to Increase the 
pixel density In the output picture may be done by Inserting the funda- 
mental discrete spectrum Into a field of zeroes and by capitalizing on a 
consequence of the similarity theorem to expand the Image. This theorem 
states that a function g(x,y) and Its transform G(fx,fy) are related by 
the expression: 


F lg(ax,by)I . |i| G 


( 10 ) 


which says tliat reducing the dimensions of the transform will Increase the 
dimensions of the Image In the function domain. This Is equivalent to 
Increasing the spacing p“l between replicates of the discrete transform, or 
decreasing the spacing between pixels. When this Is done, the true nature 
of the scene becomes apparent In Its pictorial representation. 

To illustrate the resampling operation and to display the visual 
Illusion in two dimensions. Figures 20-21 liave been created. Figure 20 
consists of rows of sinusoids extending from 0 cycles/pixel (epp) at the 
top to 0.5 epp (Nyqulst frequency) in the center to 1.0 epp at the bottom. 
Notice that the bottom half Is the mirror Image of the top half and that 
th-'s Is due to the aliasing of the entire bottom half of the Image (the 
bottom half of the "scene" bears little resemblance to its Image). 

Although the aliasing phenomenon is not the subject of this section. It is 
instructive to observe it In an Image. Interference artifacts due to 
sampling of waves having frequencies greater than 0.25 epp can be observed 
in the form of hyperbolae In the center and the four sides of the figure. 
These structures do not actually exist in the target even though they are 
readily observed in the Image. An expansion of the image performed by 
simply replicating pixels produces no changes because the high frequencies 
are retained. 

In Figure 21 the Image was expanded by a factor of two by separating 
replicates In the discrete trAnsfoinn as described above. The image con- 
tains no waves with frequencies gr^^ater than 0.25 epp, and the first-order 
interference artifacts have vani^ed. Furthermore, the image quality has 
not been degraded. 

As a rule, to Insure that digital Images will be properly visually 
Interpreted, it is necessary that the band-limited scene be sampled at a 
rate more than four or five times the band limit frequency. If this condi- 
tion is not met during data acquisition, it may be achieved by properly 
resampling the image. 


/■ 
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Some Interference artifacts will remain until the pixel spacing is 
reduced to zero, which Is the case for a continuous Image. However, the 
strongest of these are removed by an Image expansion of two times. Even 
this may be infeasible with the large images produced by current sensors 
because Fourier transforms must be computed. A compromise woold seem to 
be Interpolation using a cubic spline (Reference 15) or other tunction of 
finite size but performance superior to bilinear Interpolation. 




Figure 19. Sampled sinusoids showing Interference effects of the 

sampling function. Samples have been connected by line 
segments for clarity. From top to bottom, wavelengths 
are 40, 40/6, 40/10, 40/15, and 40/19 samples. 



Figure 20. Sinusoid target showing interference patterns as hyperbolae at center and 
in central regions of the four edges. Each line in target is a sinusoid 
in the horizontal direction. Frequency changes in vertical direction 
from 0 at top line to 1 at bottom. Note mirror image caused by aliasing 
around the middle line (at the Nyquist frequency). Image size has been 
expanded by pixel replication. 
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Figure 21. Sinusoidal target showing artifact suppression after twofjic size ex- 
pansion using the similarity theorem (sine interpolation). The few 
artifacts remaining would be removed by further interpolation to 
create a more nearly continuous image. Contrast at high frequencies 
is preserved. 
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SECTION XI 


LOW DYNAMIC RANGE INDUCED ILLUSIONS 


Quantization of real and continuous scene brightness to Integer 
values is a characteristic of all digital Images. In the illustrations 
above* the quantization Interval was 2l so that a sinusoid ranging in 
amplitude fron -1 to +1 was encoded in 256 gray levels. In real scenes 
contrast may be reduced so that only a few quanta describe the range of 
amplitudes of a harmonic component of the Image spectrum. This is commonly 
true for high frequency waves in unflltered Images, but may be true for 
any ccmiponent in a filtered image. A problem can occur if the dynamic 
range of the smnpled data is sufficiently narrow that integer truncation 
seriously modifies the true intensity values in the scene. In this event 
the quantized sinusoid is represented by a series of stacked square func- 
tions, or after sampling by the sum of truncated sequences of 6 functions. 
For Instance, a sine wave of frequency 0.125 cpp which is quantized in two 
levels may be represented by tjie sequence (...-1, -1, -1, 0, 1, 1, 1, 0, 

-1, -1, -1, 0...). The Fourier transform of this sequence is not simply 
6(s + 0.125) + 6(s - 0.125), the transform of the sinusoid, but is instead 

the sum of the transforms of the three elementary sequences (0, 0, 0, -1, 0, 

1, 0, 0...), (0, 0, -1, 0, 0, 0, 1, 0...), and (0, -1, 0, 0, 0, 0, 0, 1...), 

all of which are periodic with a wavelength of 8 samples and which add to 

produce the original sequence. The transfom of these Impulse pairs are 
themselves sinusoids having frequencies of 0,125, 0.25, and 0.375 cpp and 
the spectrum of the quantized sinusoid is 


cu 


S sin 
n*l 


(2^rn 0.125s) 


( 11 ) 


In this example, if the spectrum were quantized at the same level as 
the image, each sinusoid in the spectrum would itself be a series of alter- 
nating impulses 


J ■» 

E E- 

n=— oo 


(S(s - 0.125 mn ) 


( 12 ) 


Thus harmonic overtones are created by the act of quantization, and these 
overtones may have frequencies exceeding the Nyqulst frequency. In this 
case, the overtones will be aliased. 

Figure 22 illustrates this phenomenon on a synthetic sinusoid image. 
The figure consists of four rows, each showing a sinusoid at the left, the 
log amplitude of the Fourier transform of this sinusoid to its right, and 
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a high-pass filter applied to the sinusoid at far right (after contrast 
enhancement). The difference between the four rows is in the dynamic range 
of the sinusoid. The frequency of the wave was controlled by populating a 
single element in an otherwise empty spectrum and inverse transforming into 
the function domain. The sinusoid thus created was then scaled to the 
different amplitude ranges shown in the left column of Figure 22, The 
Fourier transforms of the quantized sinusoids are reproduced in the center 
column and clearly show the presence of harmonic overtones. 

Filtering of an image affected by. this phenomenon may modify the 
amplitudes at the fundamental and overtone frequencies differently. To 
Illustrate this, frequencies less than or equal to the fundamental fre- 
quency of the quantized sinusoid in Figure 22 were suppressed. The fil- 
tered Images, shown in the right column of Figure 22, clearly show the 
presence of energy at overtone frequencies. In images of real scenes, 
this phenomenon must result most noticeably In an Increase of high- 
frequency noise. However, all frequencies will be affected in an image 
dependent fashion, producing artifacts which are only indirectly related to 
the continuous scene. 

Quantization artifacts canno*' be removed from an Image by filtering, 
except by also removing all fundamental energy at the same frequency. 
Resampling Is likewise of no use. They are best minimized by selecting 
appropriate signal amplification prior to digital encoding. However, in 
any digital Image some quantization artifacts will be found, and the 
analyst must anticipate encountering them during careful Inspection of 
digital Images, 


Figure 22. Quantization artifacts in sinusoids. Left column shows 
sinusoidal targets. Center column shows their spectra. 
Although the target Is a simple sine wave, overtones are 
present In the transforms. Right column shows filtered 
target with power at and below the fundamental frequency 
of the sine wave suppressed. In row 1 the sinusoid has 
a brightness range of 2^ gray levels; in row 2 the 
target occupies 7 ^ gray levels; In row 3 it occupies 2^ 
gray levels; and In row 4 the brightness range is 2®. As 
the target contrast Increases more overtones appear but 
the power at the fundamental frequency Increases, For a 
continuous target an Infinite number of overtones of 
zero amplitude would be present. Contrast in filtered 
images has been adjusted for display. 
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SECTIOK XII 


SUMMARY AND CONCLUSIONS 


Transformations of Images Into other coordinate systems provide two 
advantages: more Information can be obtained from the image by simplify- 

ing the manner In which It Is displayed. Algorithms designed to operate 
In one coordinate system can be Induced to operate in virtually any 
geometry. 

Simple operations such as division provide physical Information on 
the M82 Ha emission. 

The Gull maximum entropy algorithm can be induced to provide greater 
resolutiot Improvement without a visually disturbing Increase In artifact. 

Recovery of the amplitude from the phase of the Fourier transform of 
a maximum entropy restored Image produces uncertain results. The process 
can be demonstrated on a noiseless target but may not be stable for prac- 
tical Imagery without the aid of other devices. 

A technique was presented to select the most unique me ‘phologlcal 
types from a sot of all observed types of galaxies. 

It Is possible to obtain approximate magnitudes automatically with 
a simple above-threshold search algorithm when variable nebulosity Is 
present across the scene. This Is accomplished by subtracting from the 
Image a model of the local nebulosity produced with a median filter. 

Relative color and polarization maps of R-Monocerotls Indicate a 
dusty, hollow cone shaped reflection nebula Illuminated by the star at Its 
vertex. 

Solid-state electronic cameras such as CCD’s are prevented from 
sampling in accordance with the sampling theorem by their physical design. 
The principal consequence of this Is undersampling and aliasing resulting 
In visual disruption of pictures of simple geometric targets or highly 
structured scenes. However, Image filtering may reveal these effects even 
In Images of ordinary scenes. Critical sampling or oversampling may be 
achieved by degrading the optics to produce a broad polntspread function, 
which band-limits the image projected onto the sensors. 

Properly sampled scenes may be reconstructed by interpolation of the 
discre* images. This costly process is widely overlooked, resulting in 
the Introduction of interference patterns in images. Again, simple geo- 
metric targets are most obviously affected, but filtered Images of most 
scenes are adequately displayed by conventional means. 

Quantization of Image data necessarily introduces energy at harmonic 
overtones of the image spectrum. This energy is aliased if the frequency 
of the overtones exceeds 0.5 cycles/pixel. It cannot be selectively 
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removed from the Imoge by rilterlng* and the beat my to auppreaa it is to 
maximize the amplification of the aenaor before digital encoding. The 
effect la reduced aa the dynamic range la increased. 
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